###### Appendix: Figure C2
###### Robustness to Alternative Dependent Variables
gc()
rm(list = ls())
set.seed(12345)
setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) # Note: if you are not using R Studio this command will not work, set WD to source file location manually

source("functions.R")
source("data/working_eba_function_current.R") 
require(pacman)
pacman::p_load(plm, sandwich, clubSandwich, lmtest, 
               ggplot2, numbers, dplyr, Hmisc, formattable, 
               htmltools, webshot, ggpubr, panelView, fixest, plyr, 
               scales, DataCombine)
require(webshot)
require(stringr)
date <- paste(data.table::tstrsplit(Sys.Date(), "-")[c(2:3, 1)], collapse = "")
load("data/cleaned/vat_panel_cleaned15Nov.RData")

################################################################################
### Appendix Table C2 - Alternate Dependent Variables
################################################################################

candVars <- dfCand %>% 
    filter(var_name %in% colnames(dfMerge)) %>% 
    filter(key == 1)

candVars <- as.character(candVars$var_name)
    

propMisscand <- lapply(dfMerge[, c(candVars)], 
                       FUN = function(x){sum(is.na(x))/length(x)})

dfPropmisscand <- as.data.frame(propMisscand)

dfPropmisscandlong <- dfPropmisscand %>% 
    tidyr::gather(all_of(candVars), key = "var_name", value = "propmiss")

dfCand <- dfCand %>% 
    full_join(dfPropmisscandlong, by = "var_name") 

dfCandkeep <- dfCand %>% 
    filter(propmiss < 0.2)
## Subset out high missingness and estimate

candVarskeep <- dfCandkeep$var_name[dfCandkeep$propmiss <= 0.2]
numSims <- 4943
depVars <- c(candVarskeep)
depNames <- c("Executive Corruption", 
              "Public Sector Corruption", 
              "Democratization Index", # Note, this is VAN, only non-VDEM variable here 
              "Corrupt Exchange", 
              "Legislature Corruption", 
              "Liberal Democracy Index", 
              "Accountability Index", 
              "Horizontal Accountability" , 
              "Diagonal Accountability", 
              "Neopatrimonialism")

## Rescale to make between 0 and 100
dfAltdv <- dfMerge
for(i in 1:length(depVars)){
    dfAltdv[, depVars[i]] <- rescale(dfAltdv[, depVars[i]], to = c(0, 100))
}

candRes <- eba_fit(dep.vars = depVars, 
                   dep.names = depVars,
                   samp.num = numSims, data = dfAltdv)

## Save it

save(list = c("candRes", "depVars", "depNames", "numSims", "dfCandkeep"), 
     file = "results/xnat_otherdvs.RData")

## Load the RData object and create final plots

load("results/xnat_otherdvs.RData")

incVars <- depVars
incNames <- depNames
pList <- list()
for(i in 1:length(incVars)){
    pList[[i]] <- candRes[[incVars[i]]]$plot + xlab("") + ylab("") + ggtitle(incNames[i])
}
names(pList) <- incNames

## Order for final plot: c(3, 6, 7, 8, 9, 1, 2, 4, 5, 10), this does accountability DVs then corruption DVs

pCand <- ggarrange(pList[[3]], pList[[6]], pList[[7]], 
                   pList[[8]], pList[[9]], pList[[1]], 
                   pList[[2]], pList[[4]], pList[[5]],
                   pList[[10]], common.legend = T, legend = "bottom")

pdf(file = "figures/appendix/appendix_fig_c2.pdf", 
    width = 12, height = 8)
print(pCand)
dev.off()

